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Abstract 

Ergodic dynamical systems with absolutely continuous invariant prob- 
ability measures are implemented as random-number generators for 
Monte Carlo computation. Such chaos-based Monte Carlo compu- 
tation yields sometimes unexpected dynamical dependency behavior 
which cannot be explained by the usual statistical argument. We re- 
solve the problem of its origin of this behavior by considering the effect 
of dynamical correlation on chaotic random-number generators. Fur- 
thermore, we find that superefRcient Monte Carlo computation can 
be carried out by using chaotic dynamical systems as random-number 
generators. Here superefficiency means that the expectation value of 
the square of the error decreases to as ^ with iV successive obser- 
vations for N — > oo, whereas the conventional Monte Carlo simulation 
gives the square of the error as jj. The computation speed of the super- 
efRcient case does not depend on the dimensionality s of the problems 
and, hence, it is superior to the low-discrepancy sequences yielding the 
square of the error in the known theoretical bounds. By deriving 

a necessary and sufficient condition for the superefficiency, it is shown 
that such high-performance Monte Carlo simulations can be carried 
out only if there exists a strong correlation of chaotic dynamical vari- 
ables. Numerical calculation illustrates this dynamics dependency and 
the superefficiency of chaos Monte Carlo computations. 
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1 Introduction 



The applications of physical processes to efficient computation have been re- 
cently gaining attention. Among many kinds of natural phenomena, chaos has 
a peculiar merit such that while the implementations of chaos are rather easy in 
computers or in some physical devices, it exhibits stochastic behavior as well as 
the deterministic nature. Thus, it is of great interest to investigate what kinds 
of efficient computations can be harnessed by chaos. Since chaotic phenomena 
are described by deterministic dynamical systems, the numbers generated by 
chaotic dynamical systems can be considered as not purely random numbers 
but pseudo-random numbers. So the right question must be focused on the dis- 
tinguishable things between chaos and the conventional pseudo random-number 
generators to perform some computation. While various proposals of tests for 
good random-number generators have been already extensively studied from the 
viewpoint of general aspects of random-number generations PP, such comparison 
studies between chaos and pseudo random-number generators have rarely been 
made. One exception is the study by Phatak and Rao, which showed that the lo- 
gistic map passed some statistical tests that ideal i.i.d. (independently, identically 
distributed) random-number generators must pass |2]. In this paper, we study 
chaotic (ergodic) dynamical systems as special-purpose random-number gener- 
ators for Monte Carlo methods, which are very powerful and general methods 
with many practical applications [3J E| • 

From the beginning, the key issue of the Monte Carlo method has been the 
method of generating random numbers [S] . However, no efficient chaotic random- 
number generator which is superior to the conventional pseudo random-number 
generators in Monte Carlo simulations has yet reported. The main purpose of this 
paper is to show that superefficient Monte Carlo computation can be performed 
by utilizing dynamical correlation of chaotic dynamical systems. In this super- 
efficient case, the expectation value of the square of the error in Monte Carlo 
simulations decreases to as ^ with iV successive observations for iV — > oo, 
while the corresponding value of the conventional Monte Carlo simulations de- 
creases to as for iV — > oo. The convergence speed of superefficient Monte 
Carlo simulations is therefore greatly accelerated by making use of chaotic dy- 
namical systems as random-number generators. Since any Monte Carlo compu- 
tation can be formulated as a problem of evaluating an integral over a certain 
domain, here we consider the problem of Monte Carlo integration to elucidate 
the superefficiency of chaos-based Monte Carlo computation. We do this by ana- 
lyzing how a local dynamical characteristic (thus, non-random nature) of chaotic 
random-number generators globally affects macroscopic (statistical) observable 
errors. The Monte Carlo method is essentially based on the ergodic principle: 
When random numbers {Cj}j=i- are uniformly distributed over Q (thus ergodic 
with respect to the Lebesgue measure on Q), the following identity holds accord- 
ing to the ergodic theorem (or equivalently according to the strong law of large 
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numbers). 




A sum J2j=i f(£j) can thus approximate the integral f a f(x)dx with an error in 
the order of ^= when N is a large number. In this paper, we combine this ergodic 
principle in Monte Carlo simulations with some ideas about generating nonuni- 
form random numbers by ergodic dynamical systems [HI El E] This paper is 
organized as follows. In Section 2, our Monte Carlo algorithm based on ergodic 
(chaotic) dynamical systems is explained. In Section 3, we demonstrate that 
phenomena of dynamical dependency of chaos-based Monte Carlo simulations 
are illustrated by a set of chaotic dynamical systems with unique invariant mea- 
sures. In Section 4, we clarify the necessary and sufficient condition for achieving 
the superefficiency such that the proposed chaos-based Monte Carlo simulations 
are superior to conventional Monte Carlo Simulations. Section 5 presents nu- 
merical results of such superefncient chaos-based Monte Carlo simulations for a 
one-dimensional integration problem with integrands via Chebyshev-polynomial 
expansions. In Section 6, numerical results of superefncient chaos-based Monte 
Carlo simulations for a multi-dimensional integration problem are given. Section 
7 presents numerical results of superefncient chaos-based Monte Carlo simula- 
tions for an integration problem over the infinite support (— oo, oo). A physical 
meaning of these superefncient chaos Monte Carlo computations is given in Sec- 
tion 8. And a summary and discussion are given in Section 9. 



2 Chaos-Based Monte Carlo Algorithm 

Consider a dynamical system X n+ i = F(X n ) which is ergodic with respect to 
an invariant probability measure p(x)dx with p(x) being a continuous density 
function p : Q — > R. This means that the measure dp(x) = p(x)dx is invariant 
under time evolution and, furthermore, is absolutely continuous with respect 
to the Lebesgue measure on Q. This dynamical system can thus be seen as a 
random-number generator with a sampling measure p(x)dx on the domain Q. In 
this case, according to the Birkhoff Ergodic Theorem [HI], for any function A(x) 
satisfying 

\A(x)\dx < oo, (2) 

the time average, limjv^oo jj X^o* pffi) * s ec L ua l to the space average J n A{x)dx 
for almost every Xq with respect to the probability measure p(x)dx such that 

~Afp = (A/p) = A(x)dx, (3) 

where 

i N-l 

B= hm TtE^™ ( 4 ) 

TV— >oo 1\ f—r 
i=0 
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and 

IB) = ( B(x) -p(x)dx. (5) 
Jn 

This means that the space average represented by J n A(x)dx can be computed by 
the time-average of successively generated observables B t = B(Xj) = A(XA / p(Xi) 
by generated by this chaotic dynamical system. In the same way, we can obtain 
a Monte Carlo algorithm for multiple-dimensional integrals: 

Vis = Hm 1 £ ''iff = (A/Ps) = I A{x)dx, (6) 

where p s (x) = T\i=i p( x i), dx = l\i=i dxi, and (B(x)) = j n B(x)p s (x)dx. Here, 
we assume that each dynamics Xj^ + \ = Fj(Xj^) has the same invariant proba- 
bility measure p(x)dx. Thus, any ergodic dynamical systems with respect to an 
explicit invariant probability measures p(x)dx, in principle, can serve as random- 
number generators for Monte Carlo simulations in calculating arbitrary multiple 
integrals. However, such dynamical Monte Carlo computation cannot be carried 
out unless their invariant densities p(x) are explicitly known. In recent years, 
many chaotic dynamical systems have been derived, which have explicit invariant 
measures with continuous densities |SJ E] ■ Hence these chaotic dynamical sys- 
tems can offer their applications for use as random-number generators for Monte 
Carlo simulations. Clearly, these random-number generators using chaos have 
strict time-correlation which is not desirable for ideal random-number generators. 
Suppose iV successive observations, A; t = A(XA,pi = p(Xi), i = 1, • • • , N, of 
quantities A(x) and a density function p{x) have been stored. We consider the 
expectation value of the square of the error as 

i N-1 

^^(([-^A/Pi-iA/p)} 2 )), (7) 

where the expectation of B denoted by ((B)) means an ensemble average with 
respect to the initial conditions X with a sampling measure p(x)dx. As proved 
in Appendix A, this expectation value <r(N) is given by the two-point correlation 
functions of B ri = B(Xj) as follows: 

1 2 N ?' 

a(N) = -{(B 2 ) - (B) 2 } + - £(1 - ^){(B B 3 ) - (B) 2 }, (8) 

where B(x) = A(x)/p(x). From Eq. (jHJ), we can see that cr(N) is composed of 
the statistical variance term 

v.W = j;{(B 2 ) - (B) 2 }, (9) 

which purely depends on the form of the integrand B and the dynamical corre- 
lation term 

2 N ? 
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which depends on the chaotic dynamical systems X n+ i = F(X n ) utilized as 
random-number generators. From the above formula, when there exists a nega- 
tive covariance (correlation) term such that 

N 

£{<5o^> - <5> 2 } < 0, (11) 

3=1 

the total variance can be reduced. So such negative correlated variables (anti- 
thetic variables) have been used as a variance reduction technique in the Monte 
Carlo method JT]. In the next section, we consider such a dynamical correlation 
term of chaotic random-number generators in detail. 



3 Dynamical Dependency 

To test the dynamical correlation effect in chaotic random-number generators, 
we investigate chaos-based Monte Carlo computations of simple integration prob- 
lems on the unit interval Q = [0, 1] by using chaotic dynamical systems with the 
same statistics described by a unique invariant measure. Chaotic dynamical sys- 
tems utilized as random- number generators here are listed in Table 1. In all 
of the numerical computations in this paper, we use M(= 1000) different ini- 
tial conditions X (j) for j = 1, • • • , M. We measure the numerical values of an 
empirical average of the square of the error defined by 

1 M I N-l 

V(N) = £ Mj)/PiU) - (A/P)] 2 - (12) 

3=1 i=0 

Since the numerical sampling does not obey the probability measure p(x)dx 
exactly, the observed value V(N) is not equal to the exact ensemble average 
of the square of the error cr(N) but fluctuates around cr(N). However, we can 
expect that 

V(N) -> a(N) M -f oo (13) 
because the central limit theorem holds for random variables 

1 N-l 

v 3 = [„ E Mj)/PiU) - (A/P)]\ J = 1, • • ■ , M (14) 

for M — > oo. Thus, the simulation results here does not depend on the pre- 
cise sampling information of the initial data {Xo(j)}j = \ i ... t M- Figure 1 shows 
that although Ulam-von Neumann map, Cubic map and Quintic map have the 
same invariant measure p(x)dx = — , dx , the resulting error variances of 

chaos-based Monte Carlo integration of I = Jq x 3 dx = \ exhibit discrepancy 
between them. Similar results hold for the generalized Ulam-von Neumann map 
and generalized Cubic map, both of which have the same invariant measure 
p(x)dx = — -, = s — . Figure 2 shows that there exists a significant 

' K ' *^x{l-x)(l-\x){l-±x) 
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dynamical dependency of the expectation value of the square of the errors for 
chaos-based Monte Carlo integration of / = Jq x(l — x)dx = | by using the 
same family of chaotic dynamical systems of Fig.l. Furthermore, these figures 
show that the chaotic random-number generator resulting in the most efficient 
computation depends on the integrand(the generalized Cubic map for Fig. 1 
and Ulam-von Neumann map for Fig. 2). Therefore, it is clear that these nu- 
merical results of chaotic random-number generators cannot be explained by the 
usual statistical argument such as the importance sampling. Hence, we must 
consider the dynamical effect of these chaotic random-number generators to ex- 
plain these dynamical dependency phenomena. In the next section, we consider 
an extreme case of dynamical dependency resulting in superefficient chaos-based 
Monte Carlo simulations. 



4 Condition for Superefficiency 

Here, we derive the condition for attaining the superefficiency of chaos-based 
Monte Carlo simulations such that the expectation value of the square of the error 
decreases to as ™ for N — > oo. As mentioned in Section 2, if random variables 
generated by chaotic random-number generators have a negative correlation, the 
resulting expectation value cr(N) of the square of the error can be reduced as the 
usual variance reduction technique. Furthermore, we have a stronger statement 
in the case that a chaotic dynamical system has a mixing (thus, ergodic) property 
as follows: 



Theorem 1 (Superefficiency Condition) The expectation value of the square 
of the error cr(A^) decreases to zero as as 

2 N 1 
a(N) = -— J2j{{B B j )-(Bf} = O(— )-0 for N - oo, (15) 

if and only if the relation 

oo 

<B 2 > - {Bf + 2 Y,i(B B s ) - <B> 2 } = (16) 

3=1 

is satisfied. 

The proof of the theorem is given as follows. Since cr(N) has the following identity 
from formula (JBJ): 

°( N ) = t t }, 



N N 2 . , 

3=1 



(17) 
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we have the relation 

1 oo 

a(N) = 0(—) for N ^ oo ^ (B 2 ) - (B) 2 + 2J2{(B B 3 ) - (B) 2 } = 0. 

iV 3=1 

(18) 

According to the mixing property of the chaotic dynamical systems with Lya- 
punov exponents A(> 0), we have 

(BqBj) - (B) 2 = 0[e^ x ] -> 0, for j - oo, (19) 

which assures the convergence of the infinite sums YfjLi{{BvB 3 ) — (B) 2 } and 
YfjLii{(BQBj) — (B) 2 } in the superefficiency condition (fT6j) . It is clear that 
Theorem 1 gives us a necessary and sufficient condition for the superefficiency 
of the chaos-based Monte Carlo simulations. Furthermore, Theorem 1 about the 
condition for superefficiency does not depend on the dimensionality s of integra- 
tion domains. Thus, we call the condition given in Eq. (jl6J) the superefficiency 
condition for general s-dimensional Monte Carlo integration problems. 

The superefficiency here is reminiscent of low- discrepancy sequences (or quasi- 
random numbers) generated by a deterministic algorithm, where the theoretical 

error bounds are known to be of size Qf foffi ] for s-dimensional integration prob- 
lems, which are also superior to the conventional Monte Carlo simulations with 
the error of size 0(-^=)[12j. Consequently, the theoretical values of the square of 

the error in the low-discrepancy sequences are in the order of Of ^ 1 ^? ]. Thus, 
a particular superefficient Monte Caro simulation using chaos is superior to low- 
discrepancy sequences in the sense that the expected value cr(N) of the square 
of the error in the superefficient Monte Carlo simulation is O(-^a) where the log- 
arithmic factor (lniV) 2s , which is not negligible for a large dimensionality s, is 
removed from the corresponding order of the low-discrepancy sequences. Table 
2 compares the data about the order of the square of the errors. Besides the 
differences in the order, the following differences between the two deterministic 
methods of Monte Carlo computations exist. While it is difficult for explicitly es- 
timating the error for low-discrepancy sequences, this superefficient Monte Carlo 
computations can give an exact estimate of the mean error variance as in Eq. 
(fL5|) . However, the present superefficient chaos-based Monte Carlo computations, 
which are superior to low-discrepancy sequences in higher dimension s, can be 
carried out only if there exists a negative correlation of chaotic dynamical vari- 
ables satisfying superefficiency condition ()16|) . In the following sections, we will 
present concrete examples yielding superefficient chaos Monte Carlo computa- 
tions. 

5 Superefficiency: A Case of Chebyshev Maps 

Here, we achieve superefficiency o~(N) = O(^) for a certain family of one- 
dimensional integration problems by using a set of specific chaotic dynamical 
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systems X n+ i = T p (X n ), where T P (X) is the p-th order Chebyshev polynomial 
defined by T p [cos(6)] = cos(p9) at p > 2. Examples of Chebyshev polynomials 
are given by 

T (X) = 1, T X {X) = X, T 2 (X) = 2X 2 - 1, T 3 (X) = 4X 3 - 3X, ■ • • . (20) 

It is known [7j that these Chebyshev maps T p are mixing (thus, ergodic) with 
respect to the invariant measure n J^_ x2 on t ne domain ft — [— 1, 1] for p > 2 
and they have Lyapunov exponents Inp. Note that the chaotic dynamical system 
X n+ i = T 2 (X n ) = 2X 2 — 1 is equivalent to the well-known logistic dynamical 
system Y n+ i = 41^(1 — Y n ) with an invariant measure — , dy on the unit 

interval [0, 1][6J. In the same way, chaotic dynamics caused by the third-order 
Chebyshev map X n+ i = T 3 (X n ) = 4X^ — 3X n is equivalent to the cubic chaos 
Y n+ i = Y n (3 — 4Y n ) 2 with Lyapunov exponent ln3. It is also known that a system 
of Chebyshev polynomials constitutes a complete orthonormal system satisfying 
the relations 

J l T t {x)T ] {x)p{x)dx = 5 j 1+ 2 6l '°\ (21) 
where 5jj stands for the Kronecker delta function such that 

= (22) 

Consider a one-dimensional continuous function B(x) as an integrand on the 
domain Q = [—1, 1]. By Weierstrass's approximation theorem, for an arbitrary 
e > 0, there exists a polynomial function P(x) satisfying 

\B(x) - P(x)\ < e for - 1 < x < 1. 

We can thus uniquely expand an arbitrary integrand (continuous function) B(x) 
in terms of Chebyshev polynomials (orthogonal-basis functions): 

oo oo 2 

B(x) = ]T b k T k (x) = Y -——(BT k )T k (x), 

where the coefficients are given by the formula 

2 r 1 



I B{x)T k p{x)dx = ——(BT k ). 



(1 + 6k,o) J-i 1 + 4,o 

The completeness of the system of Chebyshev orthogonal functions assures that 

b k -> for k -> oo. (23) 

Thus, it can be said that any continuous function can be well approximated in 
terms of finite numbers of Chebyshev polynomials. So we consider the Monte 
Carlo integration problem of an integrand 

L 

B(x) = J2bkT k (x). (24) 

k=0 
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Since 

= E ^T fc [T, (x)] = J2 hT kj (x), (25) 

fc=0 k=0 

we can calculate two-point correlation functions (B Bj) (= (£> (x)£> [T p j (x)])) com- 
puted by the p-th Chebyshev map T p as follows: 

L L 

(B Bj) = (b + J2 b kT k (x),b + J2 b kT k . p j(x)) 

k=l k=l 



bl + \Y} k L lf ] h p] b k j<[log p L] 
bl j>[log p L] + l 



(26) 



where [x] is the largest non-negative integer part of x. It can be seen from Eq. 
f|26|) that when p > L, the p-th Chebyshev map causes no dynamical effect such 
that <j(N) = a s (N). Hence from now on, we consider the case p < L. The 
statistical variance term is given by 

v.W = ^{{B 2 ) - (B) 2 } = i(^J> (27) 
and the dynamical correlation term of the variance of error is given by 

^(N) = -Y,(l-i r ){(BoB ] )-(B) 2 } = - Y: E (28) 

iV j=l iV iV j=l iV m=l 

The condition for the superefficiency in Eq. ()16|) is thus given in terms of the 
coefficients of the Chebyshev expansions of an integrand B(x) = Y^m=o^mT m {x) 
as follows: 

£{W?;> - (S) 2 } = E E b mp ,b m = --(E «4) < 0, (29) 

j=l j=l m=l m=l 

where the Monte Carlo simulations are computed by the p-th Chebyshev map 
T p . Therefore, for any integrand given by Chebyshev expansions as Eq. (|24jl . a 
dynamical correlation term represented by the L.H.S. of Eq. ()29|) must thus be 
negative to satisfy the superefficiency condition As is clearly seen in Eq. (J29*j) . the 
constant term 6 is irrelevant for the superefficiency Here we can give a family 
of examples that meet the superefficiency condition as follows. 

Consider the case that a normalized integrand B{x) = -^l has the following 
form 

l A(x) 
B(x) = a c + a m T p m(x) = — -y, (30) 

m=0 P\ X > 

where p > 2. Note that 

(B)=a c and (B 2 ) = a\ + \ E (31) 

/ m=0 



In this case, the following two-point correlation functions are simply given by 



l l ^ L 

(B Bi) = (a c + a m T p m(x), a c + ^ a m T pm +i(x)} = a 2 c + - ^ a m a m - h (32) 

m=0 m=0 m=l 

Thus, with the use of formula (jSJ), the expectation value of the square of the 
error <j(N) is explicitly given by the following formula: 

i L j L L 

<?( N ) = a mf - T^EE lOmOm-l- (33) 

ZJV m =0 JV !=lm=l 

Therefore, superefhciency condition (fTHjl is equivalent to the simple condition 

L 

E a ™ = 0- (34) 

m=0 

In other words, when condition (J3*4"|) is satisfied, the 0(1/N) term of c(iV) can 
be eliminated. Hence, if an integrand has a form of Eq. (j3*Uj) with the condition 
J2m=o a m — 0, the supereflicient Monte Carlo computation can be carried out by 
the p-th order Chebyshev chaos map such that the resulting expectation value 
of the square of the error are given as follows: 

*W = -^EE la m a m ^ = 0(±). (35) 

1=1 m=l 

On the other hand, even if an integrand satisfies the above condition, other 
types of chaotic dynamics, such as, Xj + i = T p /(Xj) for p' ^ p then cause the 
time-correlation functions to exhibit no dynamical effect: 

L L 

(B Bi) = (a c + E a m T p m(x),a c + E a m T p m. p n(x)) = a 2 c = (B) 2 . (36) 

771=0 771=0 

Therefore, <j(N) in this case does not so quickly decrease to for iV —>■ oo as the 
superefficient cases but normally converges to 0: 

a(N) = a s (N) = -L( ^ a 2 J = O(i). (37) 

777=0 

This means that specific chaotic (Chebyshev) dynamical systems do exist in all of 
the Chebyshev chaotic dynamical systems, which make Monte Carlo computation 
superefficient and this selection of the Chebyshev chaotic dynamical systems for 
the superefliciency depends on a form of integrand, irrespective of the fact that all 
of the Chebyshev chaotic dynamical systems has the same invariant probability 
measure (the same statistics). It can thus be said that the superefliciency of 
chaos-based Monte Carlo simulations is caused not by a statistical effect but by 
a purely dynamical effect (thus, a non-random effect) with a strong correlation 
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of chaotic variables in random-number generators. To numerically confirm the 
interplay between an integrand and chaotic dynamical systems chosen as random- 
number generators, we consider a simple integrand on Q = [—1,1] as follows: 

T _ 9r 2 4- 1 

A(x) = + 2 = [Ux) - T 2 {x)]p{x) = B(x)p(x). (38) 

Here, this integrand corresponds to the case that do — 1; Oi\ — — 1, L = 1, and 
p = 2 in Eq. (|50|l . which clearly satisfies the above condition of superefficiency. 

Thus, the random numbers generated by the second-order Chebyshev maps 
X n+ i = T 2 (X n ) are predicted to yield the expectation value of the square of the 
error, 

V T2 (N) « a T2 (N) = i (39) 

by the formula (|HHJ), while the other types of the p-th Chebyshev maps at p > 3 
are predicted to give the expectation value of square of error in the order of 

V Tp (N) « a Tp (N) = 1 for p > 2. (40) 

Figure 2 shows that numerical results coincide with our prediction about the 
occurrence of the superefficiency of chaos-based Monte Carlo simulations. The 
superefficiency achieved here remarkably contrasts with the conventional Monte 
Carlo simulations with V(N) ~ &(N) = O(jr) computed by the other p-th 
Chebyshev maps p > 2. If p in Eq. (j3*Uj) is a composite number such as p = p' k for 
an integer k(> 2), the p'-th order Chebyshev map T p i also gives a superefficiency. 
To check this numerically, we consider an integrand 

(-8a; 4 + 8x 2 + x - 1) 
A{x) = [Tx(x) - T 4 {x)}p{x) = . 

7TV 1 — 

In this case, both the second-order Chebyshev map T 2 and the fourth-order 
Chebyshev map T 4 are predicted to give superefncient results such as 

V T2 (N) ^ a T2 (N) = A (41) 

and 

y T4 (iV)^ ( x T4 (iV) = A (42) 

respectively, while the third-order Chebyshev map T 3 and the fifth-order Cheby- 
shev map T5 are predicted to give the normal behavior of the expectation value 
of the square of the error 

Vr p (N)na Tp {N) = ± for p = 3, 5. (43) 

Figure 3 confirms that numerical data coincide with our prediction about the 
selection of chaotic dynamical systems for performing super efficient chaos-based 
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Monte Carlo computations. The sampling dependency of superefficiency were 
also numerically tested. Figure 4(a) shows the numerically obtained error vari- 
ance V(N) coincides well with the exact mean value of the square of the er- 
ror cr(N) in Eq. (|H5jl for the superefficient Monte Carlo computation of an 
B (x) = T\ — + To with the sampling measure well approximated by the invari- 
ant measure p(x). Figure 4(b) shows that while there exists discrepancy between 
the numerical value and the exact value of the mean square of the error for the 
same problem with uniform sampling measure of initial data, the same kinds 
of superefficient Monte Carlo simulations are carried out. Thus, superefficient 
chaos-based Monte Carlo computation is very robust under our choice of the 
initial conditions of chaotic dynamical systems. 



6 Superefficiency: A Case of Mult i- Dimensional 
Integration 

As for the case of one-dimensional Monte Carlo integration problems, we try 
to achieve the same superefficiency in the multi-dimensional integration prob- 
lem. Let us consider the domain of integration as the s-dimensional cubic Q s = 
[— 1, 1} S . Even in this multi-dimensional case, there exists the multi-dimensional 
version of Weierstrass's approximation theorem ^Hl which guarantees that for an 
arbitrary e, there exists a polynomial function P(x) = P(x\, • • • , x s ) such that 

\B(x) - P(x)\ < e for 16 Q s , (44) 

where B(x) = B( continuous function over the domain From 

the complete orthonormal property of the Chebyshev polynomials {Ti(x)}, we 
can construct Tj 1 (xi) ■ ■ -Tj s (x s ) as an element the complete orthonormal basis 
functions spanning the functional space £ 2 (fl s ,p s ) such that 

C 2 (Q s ,p s ) = {B(x x , ■ ■ ■ ,x s )\ / \B(x x ,---,x s )\ 2 p s (x)dx 1 ---dx s < 00}. (45) 

Any continuous function B(xx, ■ ■ ■ , x n ) G £ 2 ({l s ,p s ) can thus be uniquely ex- 
panded in terms of the products of the Chebyshev polynomials as follows: 

B(x ir --,x a )= b 3i,-j. T 3iM'-- T 3.( x ')- ( 46 ) 

h,-,js>o 

Therefore, as is the one-dimensional case, superefficiency condition are given for 
any continuous function B(xi, • • • , x n ) G £ 2 (i7 s , p s ) in terms of coefficients of the 
Chebyshev expansions. Here, to present mult i- dimensional superefficiency, we 
consider a more specific family of integrands with the form, 

^1 = B(x) = B(x 1} ■ ■ ■ ,x a ) = a c + a m T pT (x l ) ■ ■ -T pT {x s ), (47) 

Ps{ x ) m =0 
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where each dynamical variable Xij associated with Xi is computed by the Pi(> 
2)-th order Chebyshev map as Xij+i = T Pi (X it j). Note that a polynomial 
T(xi, ■ ■ ■ ,x n ) = T p mi (xi) ■ ■ ■ T p m s (x s ) for p^ > 0,mj > is also an element of 
orthonormal basis-functions which constitutes a complete orthonormal system in 
the s-dimensional functional space £ 2 (f2 s , p s ). As in the one-dimensional case, by 
using the orthonormal property of the Chebyshev polynomials, we can compute 
the following two-point correlation functions: 

1 L 

{B^Bi) =a 2 c + — J2 a m a m-i- (48) 

m=l 

As will be proven in Appendix B, the expectation value of the square of the error 
(t(N) is thus given by 

I l L L 

a(A = ^ly ( E a mf - E E lOmOm-l- ( 4 9) 

m=0 1=1 m=l 

Therefore, if an integrand satisfies the superefficiency condition J2m=o a m — 0, 
which is equivalent to the relation |T6|) . s chaotic dynamical systems A ij+ i = 
T Pi (Xij) for i = 1, • • • , s again yield a superefficient Monte Carlo computation 
with the expectation value of the square of error in the order of -^2 as follows: 

°W = ~^-2 E E l*ma m -i = O(^). (50) 

1=1 m=l 

Thus, the superefficiency of chaos-based Monte Carlo computation can be 
carried out even for multi-dimensional integration problems. 

To numerically check the mult i- dimensional superefficiency of chaos-based 
Monte Carlo simulations, we consider a two-dimensional integrand as follows: 

B(x,y) = xy- (2x 2 - l)(4y 3 - 3y) = Ti(x)Ti(y) - T 2 (x)T 3 (y). 

This integrand corresponds to a c = 0,a = l,ai = — l,pi = 2,p 2 = 3 and 
L — 1 in Eq. (|4Tj) . Thus, the superefficiency condition is satisfied. Therefore, 
if we use the chaotic dynamical systems Xj + i = T 2 (Xj) and Yj + % = T 3 (Yj) as 
random-number generators for x and y respectively, superefficient Monte Carlo 
computation is predicted to be carried out such that the numerically obtained 
expectation value of the square of the error V(N) is given by 

V T2 , n (N) « v T2 ,T 3 (N) = -L (51) 

while the other chosen sets of chaotic dynamical systems yield 

Vt p1 ,t p2 (N) « a Tpl>Tp2 (N) = i (52) 

for (pl,p2) = (2,2), (3, 2), or, (3,3). Figure 5 shows that the numerical data co- 
incide well with our theory about the multi-dimensional superefficiency of chaos- 
based Monte Carlo simulations. 
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7 Superefficiency: A Case of Non-Gaussian Chaos 



Let us consider an integration problem over the infinite support Qj = (—00, 00) 
with density function Png(%)- In this case, Weierstrass's approximation theorem 
cannot be directly applied. However, when we consider a functional space 

^{Q^Png) = {f(x)\ f \f(x)\ 2 p NG (x)dx < 00}, (53) 

the chaos-based Monte Carlo integration of an integrand f{x) G £ 2 (i?/, Png) can 
also be performed without any modification of the above algorithm. The Monte 
Carlo integration over the infinite support also has many applications such as 
the pricing of exotic options in financial markets. Furthermore, the problem of 
Monte Carlo integration over the infinite support is very important from the 
theoretical point of view of Monte Carlo computation, because we cannot di- 
rectly use the usual uniform random-number generators but must consider some 
non- uniform random numbers over Usually, a procedure of Monte Carlo 
integration over the infinite support must include a transformation from the 
uniformly distributed random numbers to random numbers with non-uniform 
distributions such as the inverse method [T%j. Here without resort to a transfor- 
mation of random variables, we directly use certain ergodic dynamical systems 
with invariant probability measures over the infinite support Qi = (— 00, 00) as 
random- numbers generators over Qj. Such ergodic dynamical systems over the 
infinite support were systematically found from the multiplication formulas of 
tan(#) and its related functions [T3]. We can thus use these chaotic dynamical 
systems for this chaos-based Monte Carlo simulations over the infinite support 
Qi. Let us briefly explain ergodic mappings over Qj. As will be shown in Ap- 
pendices C and D, by using the topological conjugacy relation with Chebyshev 
mappings Tj, we can derive an infinite number of ergodic transformations 

F l (y)=h~ 1 oT,o%), (54) 



where h is a differential onto-mapping (diffeomorphism) such that h : (—00, 00 

(— 1, 1) is given by My) = sgn(y) for a > 0. These ergodic mappings 

^/i+\y\ 2a 

{Fi} have the same invariant probability measure ^H] 



pMv) = ^i + i^)- ( 55 ) 

Furthermore, as will be shown in Appendix C, we can consider the correspond- 
ing orthonormal system of functions {Pi{y) = T\ o h(y)} satisfying the same 
orthogonal relation as the Chebyshev orthogonal polynomials as follows: 

r°° (1 -\- fi n) 

/ P i (y)P j (y)pNG(y)dy = 6 iJ { + " • (56) 

These orthogonal functions constitute a complete orthonormal system of func- 
tions in £ 2 (i7/, Png)- Hence, as in constructing multidimensional integrals in 
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terms of Chebyshev polynomials where sirperefficiency condition is attained, 
chaos-based Monte Carlo simulations for a family of integrands 

L 

B(x u ■ ■ ■ ,x 8 ) = a c + E a m P P ™{xi) • • -P p ™{x s ), (57) 

where an each dynamical variable X^j associated with Xi is computed by Xij+i = 
Fp.(Xij), are superefhcient if the condition 

L 

E «m = (58) 

m=0 

is satisfied. Table 3 lists these transformations {-F)}z=i,2,3 and their dual transfor- 
mations {F*}i = i 2,3 over Qi and their related orthogonal functions {-Pz(x)}z=i,2,3 
and {P*(ic)}i=i,2,3 which were used for numerical simulations to attain the su- 
perefficiency over the infinite support. Figure 6(a) shows the graphs of er- 
godic transformations at a = 1. Figure 6(b) shows the graphs of the density 
functions Png(x) in Eq. (|76|) over Qi at several different a. In Fig. 7, we 
show that superefficient Monte Carlo computations of one-dimensional integrands 
B(x) = P 2 (x)-P 1 (x) = f=J - 5ggf and 5(x) = P|(x) -Pf (x) = - ^ 
(see Table 3) are carried out by the corresponding second-order ergodic map- 
pings F 2 (X) and F£{X) at a = 1. In Fig. 8, we show that the superefficiency 
of chaos-based Monte Carlo computation is carried out for a two-dimensional 
integral of an integrand 

1 - x 2 ) (1 - 3y 2 )sgn(?/) sgn(x)sgn(?/) 



B(x,y) = P 2 (x)P 3 (y)-P 1 (x)P 1 (y) 



[x 



2 



+ i)(i + y 2 )l VTT^VT+u 1 



(59) 

(see also Table 3) over the infinite square ^2 = (— oo, +oo) x (— oo, oo) when we 
use two different chaotic dynamical systems X n+ i = i^P^n) and Y n+ i = F^iYn) 
at a = 1 for the corresponding variables x and y of the two-dimensional integrand 
B(x,y). Hence, it is shown that chaos-based Monte Carlo computations can 
be performed to the integration problems over the infinite support by utilizing 
ergodic transformations over Qj and superefficient Monte Carlo computation can 
also be attained if the superefficiency condition is satisfied as before. 



8 Physical Meaning of Superefficiency 

Since the ergodic dynamical systems in this type of Monte Carlo computations 
generate Markov processes jTHj, we can naturally construct stochastic processes 
obeying Brownian motion from the dynamical systems. We consider a random 
variable SB{ = Bi — (B) = B(Xi) — (B(Xi)) generated by a chaotic dynamical 
system X n+ \ = F(X n ) as a velocity of a Brownian particle at t — i. Here, the 
particle position of this discrete-time Brownian motion at t = N is given by 

N-l N-l 

r(N)= E^ = E^i -<£»■ (6°) 

8=0 1=0 
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The diffusion coefficient Dd iscre te is thus obtained by the formula 

^discrete - i^J-^P^ = = ^ ~ ^ +£{ W^)-^) 2 }- 

j=i 

(61) 

In the second equality we used Eq. ©. Here, formula (J6T|) can thus be seen as 
a discrete-time version of the formula 

((r(t) 2 )) r°° 

^continuous = = L ^<^ dt ' ^ 

for continuous-time Brownian motion {r(t)} with the particle velocity u(t) sat- 
isfying the Langevin equation 

u{t) = -7«(t) + R{t) (63) 

where i?(i) is assumed to be a Gaussian process such that 

((R(t))) = 0,((R(t)R(t'))) = 2e5(t-t') (64) 

and r(t) = fQu(f)df. Hence, superefficiency condition (116)1 is simply rewritten 

as 

^discrete = ~ o ~ + XX W?;> " ^ = °- ( 65 ) 

Note that -Ddi scre t e = does n °t necessarily mean the corresponding discrete- 
time Brownian motion is dead, such as Brownian motion at zero temperature. As 
is clearly seen in Fig. 9, the discrete-time Brownian motion which corresponds to 
the superefficiency -Discrete = is still actively fluctuating around zero without 
diffusing while the discrete-time Brownian motion corresponding to the normal 
Monte Carlo simulations, -Ddi scre te > 0, has a normal diffusion behavior such 
as (r 2 (N)) = O(N). We call such non-diffusing Brownian motion ultracoherent 
Brownian motion. This ultracoherent Brownian motion (-Ddi scre t e — 0) is quite 
unlike the zero diffusion (-D cori tinuous = 0) case °f continuous-time Brownian 
motion. Hence, the physical occurrence of the ultracoherent Brownian motion 
corresponding to the superefficient chaos Monte Carlo simulation is owing to the 
discrete nature of time in chaotic dynamical systems. Whether such ultracoher- 
ent Brownian motion is observed in real physical systems would be an interesting 
problem beyond the scope of the present paper. 

9 Summary of the main results and conclusions 

In this paper we show that Monte Carlo computation of a general multi-dimensional 
integrals are performed by various ergodic dynamical systems utilized as non- 
uniform random-number generators. Such chaos-based Monte Carlo computa- 
tions exhibit a dynamical effect of random-number generations that cannot be 
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explained by the statistical arguments such as the importance sampling. By 
evaluating the expectation value cr(N) of the square of the error in terms of 
two-point correlation functions, the dynamical dependency of chaotic random- 
number generations is shown to cause a non-negligible effect in Monte Carlo 
simulations. While an effect of chaotic random-number generations is gener- 
ally seen as a variance reduction effect caused by negative covariance, we have 
shown that a superefficient Monte Carlo computation can be carried out such 
that the expectation value of the square of the error decreases to as with 
iV successive observations for iV — > oo. Therefore, a proper choice of chaotic 
random-number generators to cause superefficiency greatly accelerates the speed 
of Monte Carlo computations. We have given the necessary and sufficient con- 
dition for achieving superefficient Monte Carlo computation. Interestingly, this 
superefficiency condition has a physical meaning such that the diffusion coef- 
ficient of the corresponding discrete-time Brownian motion constructed from 
chaotic random- number generators vanishes (D 'discrete = 0), where the Brown- 
ian particles move actively in a localized region around the expectation value 
of the Monte Carlo simulation. We numerically verified the superefficiency of 
chaos Monte Carlo computations which sharply depends on the integrand and 
our choice of chaotic dynamical systems as random-number generators under var- 
ious settings of numerical integration problems. Thus we can see that observed 
dynamical dependency resulting in variance reduction in Section 3 is a partial 
realization of an extreme case of dynamical dependency, i.e., superefficiency. 

Monte Carlo methods are widely used in various computational problems 
(such as calculating financial derivative securities, communication traffic analysis, 
optimization problems, and biological computations.) Therefore, such remark- 
ably superior results from using chaotic correlation here give us to hope that they 
can effectively be applied to such hard computational problems. Furthermore, 
it is of special interest for us to investigate an possibility that some physical (or 
biological) computations with unknown mechanism such as fast protein-folding 
computation, would utilize the present kinds of superefficient Monte Carlo com- 
putation by the strong chaotic correlation of random-number generations. 

In conclusion, the present paper shows that proper chaotic (ergodic) random- 
number generators have a distinguishable effect from the conventional pseudo 
random-number generators such that chaotic random-number generators can 
greatly accelerate the speed of convergence of the mean square error in Monte 
Carlo computation (O(jj) — > O(^)) by utilizing negative correlation. Such 
superefficient Monte Carlo computation is robust under a change of the initial 
conditions, which is very promising for the application of chaotic random-number 
generators to many other computational problems, in which the expectation value 
of the square of the error is usually considered to slowly converge to as O(j^). 
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A Derivation of Eq. (JSJ) 

In this appendix, we derive formula (JHJ). 



i N-l 

*W = «b^E^-£] 2 » 

3=0 



i N-l 1 N-l 

% E B 3?)) - 2«[^ E b s ]B)) + ((W 



j=0 j=0 

2 



i N-l o N-l 

JH E <£J> + jjS E (N - j)(B oBj ) - 2(B 2 ) + (B) 

j=0 3=1 



= ^(B 2 )+^ 2 -^ 1 (B J B )-(Br 

1 2 N 

= -TriiB 2 ) ~ (B) 2 } + — £(iV - j){{BoB 5 ) - (5) 2 }- 

i=i 

In the third equality, we used the invariance property: 

(B[F{X)\) = (B) (66) 

and the following identities 

1 N-l i iV-1 

(fcE^-]»= vE(W> = ^> = ^ (67) 

JV 3=0 JV 3=0 

In the last equality, we use the identity 

2 N 1 

W2 52(N - J)(B) 2 = (B)* - -(B)\ (68) 
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B Derivation of Eq. ( 1491) 



In this appendix, we derive formula (j49j) . 
Since the following equalities 

2 E(i-^){(^)-(5) 2 } = ^D 1 -^^ 

1=1 1=1 m=l 



I L L j L L 

7^{(J2 a m) 2 - ( J2 fl m)} - TT^TT E E a ma m -ll, 
Z m=0 m=0 JVZ 1=1 m=l 



and 



1 L 



.2 



<£ 2 > - (5) 2 = ^E« m , 

Z m=0 

hold, we have for iV > L 

1 2 ^ / 

a(iv) = -{(5 2 )-(fi) 2 } + ^E(i-^){(^)-(fi) 2 } 

~ N 2 s N 2 s 2 S -W^ L , m m_i 

(=1 m=l 

N 2 s 2 s ~ l N 2 ^ ^ m ~ l 

1=1 m=l 

C Orthonormal System of Functions Related 
to the Topological Conjugacy Relation with 
Chebyshev Maps 

In this appendix, we derive various orthonormal functions from ergodic mappings 
by using the topological conjugacy relations with Chebyshev polynomials. Let 
us consider the topological conjugacy relation with the Chebyshev maps / G 

{■^m}m=0,- 

f(v) = h- 1 ofoh(y), (69) 

where / and / are onto-mappings such that / : Qj — > Qj and / : fi — > Q, 
and h is a differentiate one-to-one mapping (diffeomorphism) as h : Qj —>■ Q. 
Since h is a diffeomorphism, h preserves the ergodicity, and the mixing property 
of /; / also has a mixing (thus, ergodic) property if / has the mixing (thus, 
ergodic) property. Recall that the p-th Chebyshev maps at p > 2 have the 
invariant measures p(x)dx = ^jj=p and they satisfy the orthogonal relation with 

respect to this measure p(x)dx. Suppose that the mapping / has the invariant 
probability measure PNG(y)dy. Then, the probability preservation relation 

p(x)\dx\ = p NG (y)\dy\ (70) 
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holds. We thus obtain the density function Png(v) for / as follows: 

PN G (y) = P [h(y)}\^\. (71) 

Here we consider a system of mappings {i 7 l};=o,---,oo where Fi is defined by the 
topological conjugacy with the Z-th order Chebyshev map 

F l (y) = h- 1 oT l oh(y). (72) 

Since the Chebyshev maps {T;} at Z > 2 have the mixing property, all of the 
corresponding maps {i*}} at > 2 also have the mixing property. We therefore 
have the following theorem. 

Theorem 2 A systems of functions {i^(y)}j=o,-,oo defined as 

P l (y)^T l [h(y)} = h- 1 [F l (y)} 
constitutes an orthonormal system of functions such that 

[ P l (y)P,(y)p NG (y)dy = S l / 1 + ko) , (73) 

where Png(u) a density function defined in Eq. |77[ ). 
The proof of the theorem is given by the following identity 

/ P i (y)P j (y)p NG (y)dy = / h o F t o h~ l (x) ■ ho Fj o h~ l (x)p(x)dx 

= J T { (x) ■T j (x)p(x)dx (74) 

(i + M 

In the last equality, we use the orthonormal property of the Chebyshev system 
in Eq. (|21jl. The topological conjugacy relation with the Chebyshev maps thus 
gives the corresponding orthonormal system {P[(y)} with the density function 
Png(v) as we h as the corresponding ergodic mapping Fi(y) with the invariant 
probability measure PNc{.y)dy- 



D Non-Gaussian Chaos Mappings on the Real 
Line 

In this appendix, we derive various ergodic transformations on the real line 
(— oo, oo) via the topological conjugacy relations with the Chebyshev mappings 
Ti(x). Let us consider a diffeomorphism h : (— oo, oo) — > (—1,1) given by 

h(y) = =sgn(y) for a > 0. (75) 

1 2a 



\y\ 



20 



The topological conjugacy relation with the Chebyshev mappings 7} gives a set 
of infinite numbers of ergodic transformations Fi(y) = h^ 1 o XJ o h(y) with the 
invariant measure PNG(y)dy, where the density Png(v)^ s given by the formula 

^)=*(.)]l^l = ^5. (76) 

Note that when a = 1, the density function (JTBj) corresponds to the Cauchy den- 
sity function Png(v) — n (i+ y 2) with infinite variance. In general, density function 
(|75jl for < a < 2 does not have a finite variance. As a result, the central limit 
theorem is not applicable here and then a superposition of random variables 
obeying this density does not converge in distribution to the Gaussian distribu- 
tion but converge to the non- Gaussian Levy's stable distribution [H]. Thus, we 
call these dynamical systems non- Gaussian chaos. Here, we consider a function 
f(x) as an integrand in the functional space C 2 (f2j, Png) such as 

C\Qi,p NG ) = {f(x)\ I \f(x)\ 2 p NG (x)dx < oo}, (77) 

J Qj 

which are spanned by the orthonormal system of functions {Pi(y)}. The func- 
tional space C 2 (f2i, png) for < a < 2 is clearly different from the func- 
tional space £ 2 (f2i,e~ x2 ) which can be spanned by the complete orthonormal 
system of the Hermite polynomials Hi(x) = (—l) l e x2 -£je~ x2 with respect to 
the Gaussian measure e~ x2 dx. According to Ref. ^H], we can consider an er- 
godic transformation Fi(y) with the density (fTB|) as the multiplication formula 
of | tan(0)| «sgn[cos(0)], that is, 

Ft(\ tan(0)|«sgn[cos(0)]) = | tan(Z0)| «sgn[cos(Z0)]. (78) 

In this way, we can easily obtain a series of ergodic transformations {i 7 J}/=2,3,... 
with non-Gaussian density function over the infinite support (— oo, oo) as follows: 



(79) 

In the same way as used in the preceding appendix, orthogonal functions {Pi(y)} 
satisfying Eq. (J73*|) with respect to the density Png(u) over the infinite support 
fli = (— oo, oo) are derived as follows: 





\y 


|2a + 1 




y\ 


2a + 1 



3|y 



2a 



P (y) = 1, P 1 (y) = h(y), P 2 (y) = ^ , - , P 3 (y) = - , ^ ^ sgn(y), 

(80) 



1 2a) 2 



\y 

On the contrary, if we choose a diffeomorphism h* : fij — > Q given by 



h*(y)= , ■sgn(y), (81) 

1 + \y\ 2a 
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the corresponding ergodic transformations F* (y) are derived from the Chebyshev 
mappings TJ(y) via the topological conjugacy relation 

h*oF l *(y) = T l oh*(y). (82) 

In this case, is equivalent to a multiplication formula of | cot(6')|QSgn[cos(6 1 )] 
as 

F,*(| cot(0)|"Sgn[cos(0)]) = | cot(Z0)|»sgn[cos(Z0)]. (83) 
We thus obtain a series of ergodic transformations {F^y)} such that 

F;{y) = \\{\y\ a ~)\'sea(\y\~)^{y) = l ^^^ i °sgn[y(H 2 "-3)], • • • . 

(84) 

Correspondingly, orthogonal functions {Pj*(y)} related to the transformations 
{F*(y)}i=o r . are derived as follows: 

P*(y) = l,P*(y) = h*(y),P*(y) = ^zl,p-(y) = ly j^ °"; % gn(y), • • • . 

\y I (t + 1 2/ 1 J 2 

(85) 

In the construction of the above ergodic mappings, there exists a duality such 
that 

F l (y)Ft(z) = l, for yz = I. (86) 

Table 3 summarizes the analytical results here for the ergodic transformations 
and the related orthogonal functions, which are necessary for producing the nu- 
merical results of the superefficient chaos Monte Carlo computations over the 
infinite support Qj in Section 7. The graphs of ergodic transformations over fli 
are depicted in Fig. 6(a). 
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Table 1: Ergodic Mappings with Lyapunov exponents A on [0,1] 



Rational Mappings on [0,1] 


Densities p(x) 


A 


P O \ f \ A / -1 \/TTl ~\ T ~\ T \ 1 Z" 1 

p >(x) = Ax(l — x) (Ulam-von INeumann Map)|o| 


irW x(l — x) 


ln2 


j v >{x) = x{6 — Ax) (Cubic MapJIIj 


1 

TT^y x(l — x) 


1~ o 

mo 


fW(x) = x(5 - 20x + 16x 2 ) 2 (Quintic Map)[7j 


1 

TTyJ x(l — x) 


ln5 


f (2)/ x 4x(l-x)(l-lx)(l-rnx) iq 


1 


ln2 


J^rrA^l l-2(l+m+lm)x 2 +8lmx 3 + (l 2 +m 2 -2lm-2l 2 m~2lm 2 +l 2 m, 2 )x 4 M 


Kyfx(l-x)(l-lx)(l-mx) 


f^(x) = [x{-3 + Ax + A(l + m)x - 6(7 + m + lm)x 2 
+ 12lmx 3 + (I 2 + m 2 - 2lm - 2l 2 m - 2lm 2 + l 2 m 2 )x 4 } 2 ]/ 

[1 - 12(1 + m + lm)x 2 + 8(7 + m + I 2 + m 2 + l 2 m + lm 2 )x 3 

+ 120lmx 3 + 6(5Z 2 + 5m 2 - 26/m - 26Z 2 m - 26/m 2 + 5/ 2 m 2 )x 4 

+24(-2/ 2 - 2m 2 - 2/ 3 - 2m 3 + Aim + 7l 2 m + 7lm 2 )x 5 

+2A(Al 3 m + Aim 3 + 7l 2 m 2 - 2l 3 m 2 - 2l 2 m 3 )x 5 

+A(Al 2 + Am 2 + Al 4 + 4m 4 + 17Z 3 + 17m 3 - 8/m)x 6 

+4(-17Z 2 m - 17/m 2 - 17/ 3 m - 17/m 3 - 8/ 4 m - 8/m 4 )x 6 

+4(4/ 2 m 4 + Al 4 m 2 - 17l 3 m 2 - 17l 2 m 3 + 17/ 3 m 3 - 54/ 2 m 2 )x 6 

+24(— / 3 — m 3 — I 4 — m 4 + l 2 m + lm 2 — l 3 m — lm 3 )x 7 

+2A(l 4 m + lm 4 + Al 2 m 2 + Al 3 m 2 + Al 2 m 3 )x 7 

+2A(l 4 m 2 + l 2 m 4 — l 3 m 3 — l 4 m 3 — l 3 m 4 )x 7 

+3(3/ 4 + 3m 4 + 4/ 3 m + 4/m 3 + 4/ 4 m + 4/m 4 - 14/ 2 m 2 )x 8 

+3(-4/ 3 m 2 - 4/ 2 m 3 - 4/ 3 m 3 - 14/ 4 m 2 - 14/ 2 m 4 )x 8 

+3(4/ 4 m 3 + 4/ 3 m 4 + 3/ 4 m 4 )x 8 

+8(— l 4 m — lm 4 + l 3 m 2 + l 2 m 3 + l 4 m 2 + l 2 m 4 )x 9 

+8(-2/ 3 m 3 + l 4 m 3 + l 3 m 4 - l 4 m 4 )x 9 } [E] 


1 

K y/x{l-x)(l-lx)(l-mx) 


ln3 
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Table 2: Convergence speed of various types of Monte Carlo 
simulations with s-dimensional integration problems 



Random-Number Generators 


Variance of Error 


Standard Arithmetical Pseudo-Random Numbers (General) 


V(N) = O(i) 


Quasi-Random Numbers (General) 


V(N) = • (IniVH 


Supereflicient Chaos Monte Carlo [under Condition (1161)] 


V(N) = O(^) 


Chaos Monte Carlo (General) 


V(N) = O(i) 
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Table 3: Orthogonal functions related to exactly solvable chaos 
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Fig. 1(a): 

The expectation values of the square of the errors for the chaos-based Monte 
Carlo computations of the integral 1(A) = Jq x 3 dx = | are plotted for various 
ergodic transformations listed in Table 1 on the unit interval [0, 1]. A yellow line 
indicates a Monte Carlo computation utilizing pseudo-random numbers with a 
density p(x) = — -, 1 constructed from the uniform random-number genera- 

tors (in Fortran 90 Library) by the inverse method. Thus, it can be seen that 
chaotic random-number generators reduce the error variance from the existence 
of dynamical correlation. 
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Fig. 1(b) 

The expectation values of the square of the errors for the chaos Monte Carlo 
computations of the integral 1(A) = Jq x(1 — x)dx = | are plotted for various 
ergodic transformations listed in Table 1 on the unit interval [0, 1]. A yellow line 
indicates a Monte Carlo computation utilizing pseudo-random numbers with a 



density p(x) 



constructed from the uniform random-number genera- 



7T-^/ — x) 

tors (in Fortran 90 Library) by the inverse method. Thus, it can be seen that 
chaotic random-number generators reduce the error variance from the existence 
of dynamical correlation. 




Fig. 2 

The log-log plot of the expectation value of the square of the error for the chaos 
Monte Carlo computations of the one-dimensional integral 1(A) = J_ x p(x){Ti(x) — 
T2(x)}dx = where the p-th order Chebyshev maps are utilized as non-uniform 
random-number generators over the interval [—1, 1] for p = 2,3, 4, 5. A superef- 
ficient Monte Carlo computation is carried out for the second order Chebyshev 
map T 2 (X). 




Fig. 3 

The log-log plot of the expectation value of the square of the error for the chaos 
Monte Carlo computations of the one-dimensional integral 1(A) = J_ x p(x){Ti(x) — 
Ti(x)}dx = where the p-th order Chebyshev maps are utilized as non-uniform 
random-number generators for p = 2, 3, 4, 5. Superefficient Monte Carlo com- 
putations are carried out for the second order Chebyshev map T 2 (X) and the 
fourth order Chebyshev map T 4 (X). 




Fig. 4(a) 

The expectation values of the square of the error for superefficient Monte Carlo 
computations of the one-dimensional integral 1(A) = j\ p(x){Ti(x) — T 4 (x) + 
l}dx = 1 are plotted. The p-th order Chebyshev maps are utilized as non-uniform 
random-number generators at p = 2, 3, 4, 5. The initial conditions M = 1000 are 

generated according to a Chebyshev map with the invariant probability measure 

„f„\ — i 




Fig. 4(b) 

The expectation value of the square of the error for superefficient Monte Carlo 
computations of the one-dimensional integral 1(A) = j\ p(x){Ti(x) — T 4 (x) + 
l}dx = 1 are plotted. The p-th order Chebyshev maps are utilized as non-uniform 
random-number generators at p = 2, 3, 4, 5. The initial conditions M = 1000 are 
distributed according to the uniform pseudo-random numbers over the domain 
fi = [-l,l]. 




Fig. 5 

The log-log plot of the expectation value of the square of the error for the chaos 

Monte Carlo computations of the two-dimensional integral 1(A) = f\ J_ 1 p(x)p(y){Tx(x)Tx(y) — 

T2(x)T 3 (y)}dxdy = in the case that the p-th order Chebyshev maps are utilized 

as non-uniform random-number generators at p = 2, 3. A superefficient Monte 

Carlo computation is carried out for the second order Chebyshev map T 2 (X) for 

the x variable and the third order Chebyshev map T 3 (Y) for the y variable. 



Non-Gaussian Chaos Mappings(alpha=1) 




Fig. 6(a) 

Non-Gaussian chaos mappings for a = 1. All of the chaotic dynamical systems 
have the Cauchy distribution as the invariant probability measure. 
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Fig. 6(b) 

Density functions of the invariant measures of non-Gaussian chaos mappings are 
plotted for a = |, 1, |, 2, 4, and 10. 




Fig. 7 

The log-log plot of the expectation value of the square of the error for the chaos 
Monte Carlo computations of the one-dimensional integrals 1(A) = ff^ {P 2 (:r) — 
Pi(x)}p NG (x)dx = and 1(A) = !^{P^(x) - P?(x)}p NG (x)dx = where 
the second-order non-Gaussian chaos mappings F 2 (x) and F£(x) at a = 1 in 
Table 3 are utilized as non-uniform random-number generators over (—00,00). 
Superefficient Monte Carlo computations are carried out for the corresponding 
second-order non-Gaussian chaos maps F 2 (X) and F 2 *(X) respectively. 
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Fig. 8 

The log-log plot of the expectation value of the square of the error for the chaos 
Monte Carlo computations of the two-dimensional integral with the integrand 
1(A) = jr oo ir oo {Pi(x)P 1 (y) - P 2 (x)P 3 (y)}p NG (x)p NG (y)dxdy = in Eq.(59) 
over the two-dimensional infinite support (— oo, oo) x (— oo, oo) where the second- 
order non-Gaussian chaos mapping F 2 and the third-order non-Gaussian chaos 
mapping F s at a = 1 are utilized as non-uniform random-number generators 
over (—00,00). A superefficient Monte Carlo computation is carried out for the 
second order non-Gaussian chaos map F 2 (X) and the third order non-Gaussian 
chaos map F$(Y). 
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Fig. 9(a) 

The Brownian motions generated by the Chebyshev maps which compute the 
integral of an integrand B(x) = T\{x) — T 4 (x) are plotted. The ordinate indicates 
an empirical mean defined by S(N) = -j= Y.¥ =x Y.f=Q l [Bi(j)- {B)\ for M = 1000. 



The initial data {X^\^ = x,-,m are generated by a Chebyshev map with an invariant 
measure p{x)dx = v J^_ x2 ■ Superefficient chaos-based Monte Carlo computations 
( the second order Chebyshev map T 2 (X) and the fourth order Chebyshev map 
T^X)), corresponds to the Brownian motions showing a very coherent behavior 
around the integral (B(x)) = 0. 




Fig. 9(b) 

The magnification of Fig. 9(a) around < N < 1000. 



